The Data.

City-level information pulled from the following sources:

LA Almanac (race, age, and income) http://www.laalmanac.com/population/po38.php http://www.laalmanac.com/employment/em12.php

LA County vaccine percent by city http://publichealth.lacounty.gov/media/coronavirus/vaccine/vaccine-dashboard.htm

LA County age adjusted incidence rates http://dashboard.publichealth.lacounty.gov/covid19_surveillance_dashboard/

Data are incidence rates (proportions), vaccination rates (proportions), and characteristics by city.

Hypotheses:

All hypothesis are on the city-level and care must be taken when interpreting the results as they describe city-level variation and not individual level variation or co-variation.

Investigation of the December 2020-January 2021 peak for COVID 19 captured by a single day for incidence rates: 2021-01-07:

How do race/ethnicity, age, and income describe variation in incidence rates?

Investigation of the August 2021 peak for COVID 19 captured by a single day for incidence rates: 2021-08-18:

How do August vaccination proportions, race/ethnicity, age, income, and incidence rates for the Dec.-Jan. peak describe variation in incidence rates during the August peak for COVID 19?

Investigation of August 2021 vaccination rates for COVID 19 captured by a single day for vaccination rates: 2021-08-18:

How do race/ethnicity, age, income, and incidence rates for the Dec.-Jan. peak describe variation in August vaccination rates?

Variable Distributions and Coding

Cities

d <- readRDS("data/analysis_data.rds") %>% 
  arrange(Date.x) %>% 
  filter(Date.x == "2021-08-18")

city.populations <- data.frame(d$city, d$population.x)
names(city.populations) <- c("City", "Population_Size")
city.populations <- city.populations[order(city.populations$Population_Size, decreasing=T), ]
kable(city.populations, format = "simple", align = 'c', caption = "Cities and Population Size")

Cities and Population Size
	City	Population_Size
121	Los Angeles	4043337
95	Santa Clarita	220424
44	Glendale	206493
62	Lancaster	161570
78	Palmdale	158968
82	Pomona	157869
108	Torrance	149268
35	East Los Angeles	125269
39	El Monte	117269
33	Downey	114263
52	Inglewood	113582
114	West Covina	108235
77	Norwalk	107622
18	Burbank	107180
25	Compton	99904
101	South Gate	98155
20	Carson	93846
97	Santa Monica	92446
48	Hawthorne	91301
119	Whittier	91216
4	Alhambra	86724
61	Lakewood	80362
15	Bellflower	77735
12	Baldwin Park	76769
68	Lynwood	72047
85	Redondo Beach	68697
11	Azusa	65963
26	Covina	65851
6	Arcadia	65735
42	Florence-Graham	64705
74	Montebello	64375
81	Pico Rivera	64284
75	Monterey Park	62262
43	Gardena	61310
51	Huntington Park	59484
104	South Whittier	59222
32	Diamond Bar	57515
80	Paramount	56023
46	Hacienda Heights	55926
88	Rosemead	55350
45	Glendora	52764
89	Rowland Heights	51022
22	Cerritos	50067
56	La Mirada	49599
5	Altadena	43620
14	Bell Gardens	43071
84	Rancho Palos Verdes	42747
73	Monrovia	42681
8	Westmont	42442
92	San Gabriel	40954
57	La Puente	40697
29	Culver City	39865
115	West Hollywood	36951
23	Claremont	36484
107	Temple City	36455
13	Bell	36332
70	Manhattan Beach	35999
58	La Verne	35322
120	Willowbrook	34913
16	Beverly Hills	34520
90	San Dimas	34516
63	Lawndale	33614
111	Walnut	30532
72	Maywood	28049
21	Castaic	27191
34	Duarte	26444
102	South Pasadena	26053
91	San Fernando	24612
28	Cudahy	24347
19	Calabasas	24323
76	East San Gabriel	24036
110	Valinda	23371
100	South El Monte	22680
64	Lennox	22542
113	West Carson	22086
105	Stevenson Ranch	20966
2	Agoura Hills	20883
67	Lomita	20729
54	La Crescenta-Montrose	19801
49	Hermosa Beach	19670
96	Santa Fe Springs	18364
7	Artesia	16795
40	El Segundo	16786
112	Walnut Park	16143
37	East Rancho Dominguez	15308
47	Hawaiian Gardens	14676
79	Palos Verdes Estates	13522
93	San Marino	13277
27	Charter Oak	13144
24	Commerce	13069
60	Lake Los Angeles	12994
69	Malibu	12961
83	Quartz Hill	12906
99	Signal Hill	11797
98	Sierra Madre	10989
116	West Puente Valley	9835
71	Marina del Rey	9411
103	South San Gabriel	8848
118	Westlake Village	8360
87	Rolling Hills Estates	8113
1	Acton	7971
59	Ladera Heights	7071
10	Avocado Heights	6775
36	East Pasadena	6403
106	Sun Village	6036
55	La Habra Heights	5455
38	East Whittier	5306
30	Del Aire	4393
3	Agua Dulce	4158
66	Littlerock	4021
9	Avalon	3869
109	Val Verde	3309
31	Desert View Highlands	2493
94	San Pasqual	2035
86	Rolling Hills	1940
50	Hidden Hills	1890
65	Leona Valley	1751
41	Elizabeth Lake	1661
53	Irwindale	1459
117	West Rancho Dominguez	1359
17	Bradbury	1069

City <- d$city

Incidence Rates

January Peak: 2021-01-07

Jan.IR <- d$adj_case_14day_rate.y
Jan.IR.cat <- relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
                      ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(Jan.IR, breaks=seq(from=0, to=round(max(Jan.IR),-3), by=200), col="red", main="January Incidence Rates for Cities", xlab="Incidence Rates")
abline(v=quantile(Jan.IR, probs=c(.25, .5, .75)), lwd=2, lty=2)

August Peak: 2021-08-18

Aug.IR <- d$adj_case_14day_rate.x

hist(Aug.IR, breaks=seq(from=0, to=round(max(Aug.IR),-2), by=100), col="red", main="August Incidence Rates for Cities", xlab="Incidence Rates")

Percent Vaccinated in Each City

Aug.Vax.percent <- d$dose1_all_c_prcent
Aug.Vax.cat <- relevel(as.factor(ifelse(Aug.Vax.percent > quantile(Aug.Vax.percent, probs=.75), "High",
                      ifelse(Aug.Vax.percent < quantile(Aug.Vax.percent, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(Aug.Vax.percent, breaks=seq(from=0, to=100, by=5), col="blue", main="Percent Vaccination Rates for Cities", xlab="Percent Vaccinated")
abline(v=quantile(Aug.Vax.percent, probs=c(.25, .5, .75)), lwd=2, lty=2)

Race/Ethnicity

Proportion of Hispanic Individuals

CityHispanic <- d$Hispanic
CityHispanic.cat <- relevel(as.factor(ifelse(d$Hispanic > quantile(d$Hispanic, probs=.75), "High",
                      ifelse(d$Hispanic < quantile(d$Hispanic, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(CityHispanic, breaks=seq(from=0, to=100, by=5), main="Proportion of Hispanic Individuals", xlab="Proportion")
abline(v=quantile(CityHispanic, probs=c(.25, .5, .75)), lwd=2, lty=2)

Proportion of Black Individuals

CityBlack <- d$Black
CityBlack.cat <- relevel(as.factor(ifelse(d$Black > 10, "High","Low")), "Low")


hist(CityBlack, breaks=seq(from=0, to=100, by=5), main="Proportion of Black Individuals", xlab="Proportion")
abline(v=10, lwd=2, lty=2)
text(10, 60, "Cutoff at 10%", pos=4)

Proportion of Asian Individuals

CityAsian <- d$Asian
CityAsian.cat <- relevel(as.factor(ifelse(d$Asian > quantile(d$Asian, probs=.75), "High",
                      ifelse(d$Asian < quantile(d$Asian, probs=.5), "Low",
                             "Med"))), ref="Low")

hist(CityAsian, breaks=seq(from=0, to=100, by=5), main="Proportion of Asian Individuals", xlab="Proportion")
abline(v=quantile(CityAsian, probs=c(.25, .5, .75)), lwd=2, lty=2)

Age

Age data is captured by the number/proportion of individuals within each city in a age category. Age categories include: “Age.Under.15”,“Age.15.17”,“Age.18.24”,“Age.25.34”,“Age.35.54”,“Age.55.64”,“Age.65.”

Three groups are created for the analysis:

Young: combined age categories for “Age.Under.15”,“Age.15.17”,“Age.18.24”
Middle: age category for “Age.35.54”
Old: combined age categories for “Age.55.64”,“Age.65.”

Age <- d[,c("Age.Under.15","Age.15.17","Age.18.24","Age.25.34","Age.35.54","Age.55.64","Age.65.")]
Age <- t(apply(Age, 1, FUN=function(v) { v/sum(v) }))

CityAgeYoung <- apply(Age[,c("Age.Under.15","Age.15.17","Age.18.24")], 1, sum)
CityAgeOld <- apply(Age[,c("Age.55.64","Age.65.")], 1, sum)
CityAgeMid <- Age[,"Age.35.54"]
CityAge <- cbind(CityAgeYoung, CityAgeOld)

CityAgeOld.cat <- relevel(as.factor(ifelse(CityAgeOld > quantile(CityAgeOld, probs=.75), "High","notHigh")), "notHigh")
CityAgeYoung.cat <- relevel(as.factor(ifelse(CityAgeYoung > quantile(CityAgeYoung, probs=.75), "High","notHigh")), "notHigh")

ftable(CityAgeOld.cat, CityAgeYoung.cat)

##                CityAgeYoung.cat notHigh High
## CityAgeOld.cat                              
## notHigh                              61   30
## High                                 30    0

Household Income

HouseholdIncome <- d$Households
HouseholdIncome.cat <- relevel(as.factor(ifelse(HouseholdIncome > quantile(HouseholdIncome, probs=.75), "High",
                      ifelse(HouseholdIncome < quantile(HouseholdIncome, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(HouseholdIncome, breaks=seq(from=0, to=round(max(HouseholdIncome),-2), by=10000), main="Household Income Within Each City", xlab="Income", col="green")
abline(v=quantile(HouseholdIncome, probs=c(.25, .5, .75)), lwd=2, lty=2)

Variables for Analysis

Distributions

CityAge.m <- cbind(CityAgeOld.cat, CityAgeYoung.cat)
CityRaceEthnicty.m <- cbind(CityHispanic.cat, CityBlack.cat, CityAsian.cat)
d.analysis <- data.frame(Aug.IR, Jan.IR.cat, Aug.Vax.cat, Aug.Vax.percent, CityRaceEthnicty.m, CityAge.m, HouseholdIncome.cat)

summarytools::view(dfSummary(d.analysis, style = 'grid',
                               max.distinct.values = 10, plain.ascii = FALSE, valid.col = FALSE, headings = FALSE), method = "render")

Variable

Stats / Values

Freqs (% of Valid)

Graph

Missing

Aug.IR [numeric]

Mean (sd) : 436.5 (178.6) min < med < max: 39 < 447 < 1086 IQR (CV) : 186 (0.4)

110 distinct values

0 (0.0%)

Jan.IR.cat [factor]

1. Med 2. High 3. Low

61	(	50.4%	)
30	(	24.8%	)
30	(	24.8%	)

0 (0.0%)

Aug.Vax.cat [factor]

1. Med 2. High 3. Low

61	(	50.4%	)
30	(	24.8%	)
30	(	24.8%	)

0 (0.0%)

Aug.Vax.percent [numeric]

Mean (sd) : 63.7 (11.2) min < med < max: 19.3 < 64.9 < 83.8 IQR (CV) : 14.6 (0.2)

104 distinct values

0 (0.0%)

CityHispanic.cat [integer]

Mean (sd) : 1.7 (0.8) min < med < max: 1 < 1 < 3 IQR (CV) : 1 (0.5)

1	:	61	(	50.4%	)
2	:	30	(	24.8%	)
3	:	30	(	24.8%	)

0 (0.0%)

CityBlack.cat [integer]

Min : 1 Mean : 1.1 Max : 2

1	:	104	(	86.0%	)
2	:	17	(	14.0%	)

0 (0.0%)

CityAsian.cat [integer]

Mean (sd) : 1.8 (0.8) min < med < max: 1 < 2 < 3 IQR (CV) : 2 (0.5)

1	:	60	(	49.6%	)
2	:	30	(	24.8%	)
3	:	31	(	25.6%	)

0 (0.0%)

CityAgeOld.cat [integer]

Min : 1 Mean : 1.2 Max : 2

1	:	91	(	75.2%	)
2	:	30	(	24.8%	)

0 (0.0%)

CityAgeYoung.cat [integer]

Min : 1 Mean : 1.2 Max : 2

1	:	91	(	75.2%	)
2	:	30	(	24.8%	)

0 (0.0%)

HouseholdIncome.cat [factor]

1. Med 2. High 3. Low

61	(	50.4%	)
30	(	24.8%	)
30	(	24.8%	)

0 (0.0%)

Generated by summarytools 0.9.9 (R version 4.1.0)
2021-11-04

Correlation

# need to make sure low is lowest when setting as.numeric (to make cor more clear)





d.num <- data.frame(Aug.IR, Jan.IR.cat, Aug.Vax.cat,
                    CityHispanic.cat, as.numeric(CityBlack.cat), as.numeric(CityAsian.cat), 
                    as.numeric(CityAgeOld.cat), as.numeric(CityAgeYoung.cat),
                    HouseholdIncome.cat) %>% 
  mutate(across(c(Jan.IR.cat, Aug.Vax.cat, CityHispanic.cat, HouseholdIncome.cat), 
                ~ as.numeric(relevel(.x, ref = "Low"))))

cormat <- cor(d.num, use="complete.obs", method="pearson")
corrplot(cormat, type="upper", order="original",
          col=brewer.pal(n=8, name="RdYlBu"),
          title = "",
          addCoef.col = "black",
          tl.cex=.5, number.cex=.5)

Pairwise Trend Plots

January IR vs. Proportion of Hispanic Individuals

CityHispanic.cont <- as.numeric(relevel(CityHispanic.cat, "Low"))
coef <- summary(lm(Jan.IR~CityHispanic.cont))$coef
ggplot(d.analysis, aes(x=CityHispanic.cont, y=Jan.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=3000, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=2500, label=paste("p-value", "=", round(coef[2,4],25))) +
  scale_x_continuous(name ="Proportion of Hispanic Individuals", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("Jan Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

January IR vs. Household Income

coef <- summary(lm(Jan.IR~as.numeric(relevel(HouseholdIncome.cat, "Low"))))$coef
ggplot(d.analysis, aes(x=as.numeric(relevel(HouseholdIncome.cat, "Low")), y=Jan.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=3000, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=2500, label=paste("p-value", "=", round(coef[2,4],23))) +
  scale_x_continuous(name ="Household Income", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("January Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

August IR vs. January IR

coef <- summary(lm(Aug.IR~Jan.IR))$coef
ggplot(d.analysis, aes(x=Jan.IR, y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=3000, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=3000, y=800, label=paste("p-value", "=", round(coef[2,4],6)))+
  xlab("January Incidence Rates")+ylab("August Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

August IR vs. Proportion of Hispanic Individuals

CityHispanic.cont <- as.numeric(relevel(CityHispanic.cat, "Low"))
coef <- summary(lm(Aug.IR~CityHispanic.cont))$coef
ggplot(d.analysis, aes(x=CityHispanic.cont, y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=800, label=paste("p-value", "=", round(coef[2,4],2))) +
  scale_x_continuous(name ="Proportion of Hispanic Individuals", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("August Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

August IR vs. Household Income

coef <- summary(lm(Aug.IR~as.numeric(relevel(HouseholdIncome.cat, "Low"))))$coef
ggplot(d.analysis, aes(x=as.numeric(relevel(HouseholdIncome.cat, "Low")), y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=800, label=paste("p-value", "=", round(coef[2,4],2))) +
  scale_x_continuous(name ="Household Income", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("August Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

August IR vs. August Percent Vaccinated

coef <- summary(lm(Aug.IR~Aug.Vax.percent))$coef
ggplot(d.analysis, aes(x=Aug.Vax.percent, y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=80, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=80, y=800, label=paste("p-value", "=", round(coef[2,4],2)))+
  xlab("August Vaccination Percent")+ylab("August Incidence Rates")

## `geom_smooth()` using formula 'y ~ x'

January Incident Rates:

Regression Model

reg.Jan.noVax <- lm(Jan.IR ~ HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)

Regression summary

kable(summary(reg.Jan.noVax)$coef, digits=3, format = "simple", align = 'c', caption = "January Regression")

January Regression
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	1087.329	431.053	2.522	0.013
HouseholdIncome.catHigh	-983.811	171.194	-5.747	0.000
HouseholdIncome.catLow	783.304	163.598	4.788	0.000
CityRaceEthnicty.mCityHispanic.cat	97.188	81.787	1.188	0.237
CityRaceEthnicty.mCityBlack.cat	65.163	165.618	0.393	0.695
CityRaceEthnicty.mCityAsian.cat	67.657	73.121	0.925	0.357
CityAge.mCityAgeOld.cat	-107.726	146.655	-0.735	0.464
CityAge.mCityAgeYoung.cat	444.903	164.853	2.699	0.008

#stargazer(reg.Jan.noVax, type = 'html', ci=TRUE, ci.level=0.95, single.row = T)

Regression ANOVA

\(R^2 =\) 0.62

kable(anova(reg.Jan.noVax), digits=3, format = "simple", align = 'c', caption = "January Regression ANOVA")

January Regression ANOVA
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
HouseholdIncome.cat	2	60973275.2	30486637.6	87.403	0.000
CityRaceEthnicty.m	3	542565.2	180855.1	0.518	0.670
CityAge.m	2	3052626.7	1526313.3	4.376	0.015
Residuals	113	39415190.7	348807.0	NA	NA

August Incident Rates:

Vaccination Only

Regression Model

reg.Aug.vax <- lm(Aug.IR ~ Aug.Vax.cat)

Regression summary

kable(summary(reg.Aug.vax)$coef, digits=3, format = "simple", align = 'c', caption = "August Regression")

August Regression
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	450.118	21.753	20.692	0.000
Aug.Vax.catHigh	-109.075	37.886	-2.879	0.005
Aug.Vax.catLow	54.182	37.886	1.430	0.155

Regression ANOVA

\(R^2 =\) 0.11

kable(anova(reg.Aug.vax), digits=3, format = "simple", align = 'c', caption = "August Regression ANOVA")

August Regression ANOVA
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Aug.Vax.cat	2	422581	211290.5	7.32	0.001
Residuals	118	3406082	28865.1	NA	NA

Full Model

Regression Model

reg.Aug <- lm(Aug.IR ~ Aug.Vax.cat + Jan.IR.cat + HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)

Regression summary

kable(summary(reg.Aug)$coef, digits=3, format = "simple", align = 'c', caption = "August Regression")

August Regression
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	298.772	114.632	2.606	0.010
Aug.Vax.catHigh	-53.373	42.644	-1.252	0.213
Aug.Vax.catLow	43.507	41.916	1.038	0.302
Jan.IR.catHigh	-37.110	47.008	-0.789	0.432
Jan.IR.catLow	-218.582	49.835	-4.386	0.000
HouseholdIncome.catHigh	59.182	48.820	1.212	0.228
HouseholdIncome.catLow	-67.747	48.738	-1.390	0.167
CityRaceEthnicty.mCityHispanic.cat	73.615	25.706	2.864	0.005
CityRaceEthnicty.mCityBlack.cat	88.676	44.028	2.014	0.046
CityRaceEthnicty.mCityAsian.cat	-5.201	20.844	-0.250	0.803
CityAge.mCityAgeOld.cat	-70.434	39.484	-1.784	0.077
CityAge.mCityAgeYoung.cat	58.694	45.717	1.284	0.202

Regression ANOVA

\(R^2 =\) 0.31

kable(anova(reg.Aug), digits=3, format = "simple", align = 'c', caption = "August Regression ANOVA")

August Regression ANOVA
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Aug.Vax.cat	2	422581.0	211290.50	8.714	0.000
Jan.IR.cat	2	233137.6	116568.79	4.808	0.010
HouseholdIncome.cat	2	98384.8	49192.40	2.029	0.136
CityRaceEthnicty.m	3	287587.3	95862.42	3.954	0.010
CityAge.m	2	144166.4	72083.21	2.973	0.055
Residuals	109	2642806.1	24245.93	NA	NA

August Vaccination Rates:

Regression Model with August Vaccination

reg.Aug.vax <- lm(Aug.Vax.percent ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)

Regression summary

kable(summary(reg.Aug.vax)$coef, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression")

August Vaccination Regression
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	59.663	6.915	8.627	0.000
Jan.IR.catHigh	-1.763	2.791	-0.632	0.529
Jan.IR.catLow	0.389	2.866	0.136	0.892
HouseholdIncome.catHigh	3.236	2.958	1.094	0.276
HouseholdIncome.catLow	-0.898	2.940	-0.305	0.761
CityRaceEthnicty.mCityHispanic.cat	3.718	1.521	2.445	0.016
CityRaceEthnicty.mCityBlack.cat	0.993	2.651	0.375	0.709
CityRaceEthnicty.mCityAsian.cat	3.690	1.176	3.137	0.002
CityAge.mCityAgeOld.cat	-2.160	2.381	-0.907	0.366
CityAge.mCityAgeYoung.cat	-6.132	2.681	-2.287	0.024

Regression ANOVA

kable(anova(reg.Aug.vax), digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression ANOVA")

August Vaccination Regression ANOVA
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Jan.IR.cat	2	1976.111	988.056	11.094	0.000
HouseholdIncome.cat	2	1059.438	529.719	5.948	0.004
CityRaceEthnicty.m	3	1696.679	565.560	6.350	0.001
CityAge.m	2	488.550	244.275	2.743	0.069
Residuals	111	9885.456	89.058	NA	NA

Multnomial Regression Model with August Vaccination

reg.Aug.vax <- lm(Aug.Vax.percent ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
reg.Aug.vax <- multinom(Aug.Vax.cat ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)

## # weights:  33 (20 variable)
## initial  value 132.932087 
## iter  10 value 78.668921
## iter  20 value 76.286634
## iter  30 value 76.182929
## iter  40 value 76.181907
## final  value 76.181903 
## converged

Regression summary for High Vaccination Rates

coef.high <- summary(reg.Aug.vax)$coefficients["High",]
OR.high <- exp(coef.high)
se.high <- summary(reg.Aug.vax)$standard.errors["High",]
z.high <- coef.high/se.high
p.value.high <- (1 - pnorm(abs(z.high), 0, 1)) * 2
r.high <- data.frame(OR.high, coef.high, se.high, z.high, p.value.high)
kable(r.high, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression for High Vax Rates")

August Vaccination Regression for High Vax Rates
	OR.high	coef.high	se.high	z.high	p.value.high
(Intercept)	459743.494	13.038	1.130	11.538	0.000
Jan.IR.catHigh	0.000	-16.605	0.000	-12741324.270	0.000
Jan.IR.catLow	6.620	1.890	0.742	2.548	0.011
HouseholdIncome.catHigh	0.783	-0.245	0.834	-0.294	0.769
HouseholdIncome.catLow	2.879	1.058	1.343	0.787	0.431
CityRaceEthnicty.mCityHispanic.cat	1.589	0.463	0.413	1.122	0.262
CityRaceEthnicty.mCityBlack.cat	0.000	-15.583	1.130	-13.790	0.000
CityRaceEthnicty.mCityAsian.cat	1.028	0.028	0.392	0.071	0.944
CityAge.mCityAgeOld.cat	1.604	0.473	0.682	0.693	0.488
CityAge.mCityAgeYoung.cat	1.068	0.065	1.265	0.052	0.959

Regression summary for Low Vaccination Rates

coef.low <- summary(reg.Aug.vax)$coefficients["Low",]
OR.low <- exp(coef.low)
se.low <- summary(reg.Aug.vax)$standard.errors["Low",]
z.low <- coef.low/se.low
p.value.low <- (1 - pnorm(abs(z.low), 0, 1)) * 2
r.low <- data.frame(OR.low, coef.low, se.low, z.low, p.value.low)
kable(r.low, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression for Low Vax Rates")

August Vaccination Regression for Low Vax Rates
	OR.low	coef.low	se.low	z.low	p.value.low
(Intercept)	3.841	1.346	2.349	0.573	0.567
Jan.IR.catHigh	0.791	-0.234	0.942	-0.249	0.803
Jan.IR.catLow	4.193	1.433	1.110	1.291	0.197
HouseholdIncome.catHigh	0.533	-0.629	1.031	-0.610	0.542
HouseholdIncome.catLow	1.750	0.560	0.934	0.600	0.549
CityRaceEthnicty.mCityHispanic.cat	0.377	-0.976	0.574	-1.701	0.089
CityRaceEthnicty.mCityBlack.cat	0.861	-0.149	0.708	-0.211	0.833
CityRaceEthnicty.mCityAsian.cat	0.081	-2.516	0.799	-3.151	0.002
CityAge.mCityAgeOld.cat	2.407	0.878	0.853	1.029	0.303
CityAge.mCityAgeYoung.cat	3.720	1.314	0.807	1.629	0.103

Model params over time

Full model fit by date
{r} adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat
No vaccination lag time in this example

Vaccination rates

Modeled as a continuous variable. Plots include dates > 04/01/2021

d <- readRDS("data/analysis_data.rds") 
lac_summarytable_dash <- read.csv("data/LA_County_Covid19_CSA_14day_case_death_table_Clean.csv")

# data clean ----
dout <- d %>% 
  mutate(Date.x = ymd(Date.x), 
         Vax.percent = dose1_all_c_prcent, 
         Jan.IR =  as.numeric(adj_case_14day_rate.y),
         Jan.IR.cat =  relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
                                                ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
                                                       "Med"))), ref="Low")) %>%
  dplyr::select(Date.x, city, adj_case_14day_rate.x, 
                Vax.percent, Jan.IR.cat, 
                HouseholdIncome.cat, CityHispanic.cat, CityBlack.cat, 
                CityAsian.cat, CityAgeOld.cat, CityAgeYoung.cat)




test <- dout %>% 
  # dplyr::filter(Date.x > "2021-03-01") %>% 
  tidyr::nest(data = -Date.x) %>% 
  dplyr::mutate(fit = purrr::map(data, ~lm(adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat, data = .x)),
                tidied = purrr::map(fit, ~broom::tidy(.x), quality = purrr::map(fit, ~glance(.x)))) %>% 
  dplyr::select(-data, -fit) %>% 
  tidyr::unnest(tidied) %>% 
  # dplyr::filter(grepl("Vax", term)) %>%
  arrange(Date.x)



# Vaccination rates ---- 
# after april 2021 (vaccination rates too low prior)
p1 <- ggplot(test %>% filter(grepl("Vax.percent", term), Date.x > "2021-04-01"), aes(Date.x, estimate)) + 
  geom_point() + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
  # geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County", 
                      Date > "2021-04-01") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Race

p1 <- ggplot(test %>% filter(term %in% c("CityAsian.catHigh", "CityBlack.catHigh", "CityHispanic.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Household income

p1 <- ggplot(test %>% filter(term %in% c("HouseholdIncome.catHigh", "HouseholdIncome.catLow")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.2))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

Age category

p1 <- ggplot(test %>% filter(term %in% c("CityAgeOld.catHigh", "CityAgeYoung.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

January incidence rates

p1 <- ggplot(test %>% filter(term %in% c("Jan.IR.catHigh", "Jan.IR.catMed")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
# ggsave(p1, width = 10, height = 8, file = "~/Desktop/test1.png")
p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

Model params over time (vax lag)

Vaccination rates

Modeled as a continuous variable. Plots include dates > 04/01/2021

d <- readRDS("data/analysis_data_14vax_lag.rds") 
lac_summarytable_dash <- read.csv("data/LA_County_Covid19_CSA_14day_case_death_table_Clean.csv")

# data clean ----
dout <- d %>% 
  mutate(Date.x = ymd(Date.x), 
         Vax.percent = dose1_all_c_prcent, 
         Jan.IR =  as.numeric(adj_case_14day_rate.y),
         Jan.IR.cat =  relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
                                                ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
                                                       "Med"))), ref="Low")) %>%
  dplyr::select(Date.x, city, adj_case_14day_rate.x, 
                Vax.percent, Jan.IR.cat, 
                HouseholdIncome.cat, CityHispanic.cat, CityBlack.cat, 
                CityAsian.cat, CityAgeOld.cat, CityAgeYoung.cat)




test <- dout %>% 
  # dplyr::filter(Date.x > "2021-03-01") %>% 
  tidyr::nest(data = -Date.x) %>% 
  dplyr::mutate(fit = purrr::map(data, ~lm(adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat, data = .x)),
                tidied = purrr::map(fit, ~broom::tidy(.x), quality = purrr::map(fit, ~glance(.x)))) %>% 
  dplyr::select(-data, -fit) %>% 
  tidyr::unnest(tidied) %>% 
  # dplyr::filter(grepl("Vax", term)) %>%
  arrange(Date.x)



# Vaccination rates ---- 
# after april 2021 (vaccination rates too low prior)
p1 <- ggplot(test %>% filter(grepl("Vax.percent", term), Date.x > "2021-04-01"), aes(Date.x, estimate)) + 
  geom_point() + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
  # geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County", 
                      Date > "2021-04-01") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Race

p1 <- ggplot(test %>% filter(term %in% c("CityAsian.catHigh", "CityBlack.catHigh", "CityHispanic.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Household income

p1 <- ggplot(test %>% filter(term %in% c("HouseholdIncome.catHigh", "HouseholdIncome.catLow")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.2))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

Age category

p1 <- ggplot(test %>% filter(term %in% c("CityAgeOld.catHigh", "CityAgeYoung.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

January incidence rates

p1 <- ggplot(test %>% filter(term %in% c("Jan.IR.catHigh", "Jan.IR.catMed")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
# ggsave(p1, width = 10, height = 8, file = "~/Desktop/test1.png")
p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

## Warning: Removed 1 row(s) containing missing values (geom_path).

COVID-19 LA County City-level Analysis

David Conti

04 November 2021

The Data.

Hypotheses:

Variable Distributions and Coding

Cities

Incidence Rates

January Peak: 2021-01-07

August Peak: 2021-08-18

Percent Vaccinated in Each City

Race/Ethnicity

Proportion of Hispanic Individuals

Proportion of Black Individuals

Proportion of Asian Individuals

Age

Household Income

Variables for Analysis

Distributions

Correlation

Pairwise Trend Plots

January IR vs. Proportion of Hispanic Individuals

January IR vs. Household Income

August IR vs. January IR

August IR vs. Proportion of Hispanic Individuals

August IR vs. Household Income

August IR vs. August Percent Vaccinated

January Incident Rates:

Regression Model

Regression summary

Regression ANOVA

August Incident Rates:

Vaccination Only

Regression Model

Regression summary

Regression ANOVA

Full Model

Regression Model

Regression summary

Regression ANOVA

August Vaccination Rates:

Regression Model with August Vaccination

Regression summary

Regression ANOVA

Multnomial Regression Model with August Vaccination

Regression summary for High Vaccination Rates

Regression summary for Low Vaccination Rates

Model params over time

Vaccination rates

Race

Household income

Age category

January incidence rates

Model params over time (vax lag)

Vaccination rates

Race

Household income

Age category

January incidence rates